Alpha Diversity
2021-07-08
library(knitr)
options(max.print="75")
opts_chunk$set(echo = TRUE, cache = TRUE, warning = FALSE, prompt=FALSE, tidy = TRUE, comment = NA, message = FALSE)
opts_knit$set(width=75)
# Loading library here
library(tidyverse)
library(ggpubr)
library(MASS)
library(segmented)
library(agricolae)
library(car)
library(effectsize)
library(rstatix)
library(patchwork)
# set path
source("assets/utils.R")Initial setup
Setup color index
# Exposure level
Expo_color <- c("#1b9e77", "#d95f02", "#8da0cb")
names(Expo_color) <- c("Low", "Medium", "High")
print(Expo_color) Low Medium High
"#1b9e77" "#d95f02" "#8da0cb" #plot(1:3, 1:3, col = Expo_color, pch = 16, cex = 4)
# Ethnicity
Ethnicity_color <- gg_color_hue(3)
names(Ethnicity_color) <- c("SANEMA", "YEKWANA", "Visitors")
print(Ethnicity_color) SANEMA YEKWANA Visitors
"#F8766D" "#00BA38" "#619CFF" Sanemas <- c("Chajuranha", "Mosenahanha", "Kuyuwininha", "Shianana-Jiyakwanha", "Washudihanha", "Sudukuma")
Yekwana <- c("Kanarakuni", "Fiyakwanha")Import the data
# Alpha diversities
alpha = readRDS("../data/alpha_imported_all.RDS")# Mapping
## This mapping files include all villagers and baseline of visitors of all body sites
mt_ms = readRDS("../data/Mapping_MS_Expo.Rds")
Mapping_work <- mt_ms %>% mutate(Expo = ifelse(Ethnicity=="SANEMA", "Low", ifelse(Village=="Fiyakwanha", "Low", ifelse(SampleGroup=="Visitors", "High", "Medium"))) %>% factor(levels = c("Low", "Medium", "High")), Age_num = as.numeric(Age), Age_num = ifelse(is.na(Age_num), as.numeric(gsub(pattern = "_months", replacement = "", Age, ignore.case = T))/12, Age_num), Age_grp1 = ifelse(Age_num>=18, "Adults", "Children"))
Mapping_work$Age_num[Mapping_work$Age=="4_Days"] <- 0
Mapping_work$Age_num[Mapping_work$Age=="1_day"] <- 0Set up the Analysis
dat_alpha <- merge(alpha, Mapping_work, by = 1) %>% mutate(Age_grp2 = case_when(Age_num<=3 ~ "Age_0-3", Age_num<=8 ~ "Age_3-8", Age_num<18 ~ "Age_8-18", T ~ "Adults"))Effects of exposure level
Indices <- colnames(alpha[c(2:5)])
Yrs <- c("2015", "2016")
BodySites <- c("Feces", "Mouth", "Nose", "Right_Arm", "Right_Hand")Initial test on feces (no need to run)
We first check the faith_pd which described the overall complexity of the microbial community.
BB = "Feces"
AA = "faith_pd"
# Year 2015
YY = "2015"
Work_dat <- dat_alpha %>% dplyr::filter(Body_Site==BB, Year==YY, SampleGroup=="Villagers")
# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Expo + Gender + Age_num:Expo + Age_num:Gender + Expo:Gender, data = Work_dat)
# Obtain the null model with no predictors, this is necessary for forward selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)
# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~ Age_num + Expo + Gender + Age_num:Expo + Age_num:Gender + Expo:Gender)
anova(regforward.out)The best model from 2015 includes Expo, Age, Expo:Age interaction and Gender (in decreasing importance).
# Year 2016
YY = "2016"
Work_dat <- Work %>% dplyr::filter(Body_Site==BB, Year==YY, SampleGroup=="Villagers")
# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Expo + Gender + Age_num:Expo + Age_num:Gender + Expo:Gender, data = Work_dat)
# Obtain the null model with no predictors, this is necessary for forward selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)
# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~ Age_num + Expo + Gender + Age_num:Expo + Age_num:Gender + Expo:Gender)
anova(regforward.out)The best model from 2016 includes Expo, Age, Expo:Age interaction (in decreasing importance), but Gender was not there.
So we will fitting the data with Expo, Age, Expo:Age interaction, and Gender for both years to keep unity.
# Final model
lm_models <- list()
for (YY in Yrs){
Work_dat <- Work %>% dplyr::filter(Body_Site==BB, Year==YY, SampleGroup=="Villagers")
fit.final <- lm(as.formula(paste0(AA, " ~ Age_num + Expo + Age_num:Expo + Gender")), data = Work_dat)
print(paste0("--------Below is ANOVA table for year ", YY))
anova(fit.final) %>% print()
print(paste0("--------Below is coefficient for year ", YY))
summary(fit.final) %>% print()
lm_models[[YY]] <- fit.final
}Next we perform the broken stick regression based on the selected final model
for (YY in Yrs){
#YY = "2016"
fit.seg <- segmented(lm_models[[YY]], seg.Z = ~ Age_num, psi = c(3))
summary(fit.seg)
plot(fit.seg$fitted.values, fit.seg$residuals)
plot(fit.seg$model$faith_pd, fit.seg$fitted.values)
plot(fit.seg$model$Age_num, hatvalues(fit.seg))
plot(fit.seg$model$Age_num, cooks.distance(fit.seg))
plot(hatvalues(fit.seg), fit.seg$residuals)
}It seems from the residual plots that the point with faith_pd smaller than 10 is dragging the curve down and causing non-random residual.
Work %>% dplyr::filter(Body_Site==BB, faith_pd<12, SampleGroup=="Villagers")It does appears that many baby are outliers, so re-do the analysis without baby.
Modeling without babies
Overview
for (bb in BodySites){
for (aa in Indices){
Work_dat <- dat_alpha %>% dplyr::filter(Body_Site==bb, SampleGroup=="Villagers", Age_num>=1)
p <- ggplot(Work_dat, aes(x = Age_num, y = get(aa), color = Expo, shape = Gender)) +
geom_point(size = 1.5) +
scale_shape_manual(values = c(16, 21)) +
scale_color_manual(breaks = names(Expo_color), values = Expo_color) +
facet_grid(.~Year) +
theme_jw() + theme(aspect.ratio = 0.8, panel.background = element_rect(color = "black")) +
labs(x = "Age (Years)", y = aa, title = bb)
print(p)
}
}
#### Exploratory Statisticals ##### Children Building full multiple linear regression models with interaction terms for children
alpha_youngchild_lm_models = list()
lm.stats.tbl = list()
for (bb in BodySites){
for (aa in Indices){
for (yy in Yrs){
Work_dat <- dat_alpha %>% dplyr::filter(Body_Site==bb, Year==yy, SampleGroup=="Villagers", Age_num>=1, Age_num<18)
lm.model <- lm(paste0(aa, " ~ Age_num + Expo + Gender + Age_num:Expo + Age_num:Gender + Expo:Gender"), data = Work_dat, contrasts = list(Expo = "contr.sum", Gender = "contr.sum"))
lm.stats = Anova(lm.model, type = 3) %>% # Type 3 ss
merge(., effectsize::eta_squared(., ci = 0.95), by.x = 0, by.y = 1, all = T) %>% # effect size estimates
rename(Factors = Row.names) %>% mutate(Year = yy, Body_Site = bb, Indice = aa)
alpha_youngchild_lm_models[[bb]][[aa]][[yy]] = lm.model
lm.stats.tbl[[paste(bb, aa, yy, sep = "-")]] = lm.stats
}
}
}
lm.stats.tbls = do.call("rbind", lm.stats.tbl) %>% rename(p_value = `Pr(>F)`) %>% mutate(q_value = p.adjust(p_value, method = "fdr"))
# full model (multiple linear regression) statistics for alpha diversity in children
#write.table(lm.stats.tbls, file = "../output/tables/alpha_lm_stats_children.txt", sep = "\t", row.names = F, quote = F)
lm.stats.tbls.2 = lm.stats.tbls %>% dplyr::filter(!is.na(Eta2_partial)) %>% mutate(Factors = factor(Factors, levels = c("Age_num", "Gender", "Expo", "Age_num:Expo", "Age_num:Gender", "Expo:Gender"))) %>% arrange(Body_Site, Indice, Year, Factors)
# p values
ggplot(lm.stats.tbls.2, aes(x = Factors, y = q_value, color = Body_Site)) +
geom_point() +
geom_path(aes(group = Year, linetype = Year)) +
facet_grid(Body_Site~Indice) +
scale_y_continuous(breaks = c(0, 0.05, 0.1, 0.2, 1.0)) +
theme_jw(base_size = 12) + theme(aspect.ratio = 0.4, axis.text.x = element_text(angle = 90, hjust = 1), panel.grid.minor = element_blank()) +
ggtitle("FDR adjusted p-value for full models of alpha diversity in Children")
# eta partial
ggplot(lm.stats.tbls.2, aes(x = Factors, y = Eta2_partial, color = Body_Site)) +
geom_point() +
geom_path(aes(group = Year, linetype = Year)) +
facet_grid(Body_Site~Indice, scales = "free_y") +
#scale_y_continuous(breaks = c(0, 0.05, 0.1, 0.2, 1.0)) +
theme_jw() + theme(aspect.ratio = 0.4, axis.text.x = element_text(angle = 90, hjust = 1)) +
ggtitle("Effect size (eta partial) for full models of alpha diversity in Children")
No strong interactions obsereved among children
all ages
Building full multiple linear regression models with interaction terms
alpha_full_lm_models = list()
lm.stats.tbl = list()
for (bb in BodySites){
for (aa in Indices){
for (yy in Yrs){
Work_dat <- dat_alpha %>% dplyr::filter(Body_Site==bb, Year==yy, SampleGroup=="Villagers", Age_num>=1)
lm.model <- lm(paste0(aa, " ~ Age_num + Expo + Gender + Age_num:Expo + Age_num:Gender + Expo:Gender"), data = Work_dat, contrasts = list(Expo = "contr.sum", Gender = "contr.sum"))
lm.stats = Anova(lm.model, type = 3) %>% # Type 3 ss
merge(., effectsize::eta_squared(., ci = 0.95), by.x = 0, by.y = 1, all = T) %>% # effect size estimates
rename(Factors = Row.names) %>% mutate(Year = yy, Body_Site = bb, Indice = aa)
alpha_full_lm_models[[bb]][[aa]][[yy]] = lm.model
lm.stats.tbl[[paste(bb, aa, yy, sep = "-")]] = lm.stats
}
}
}
lm.stats.tbl = do.call("rbind", lm.stats.tbl)
#o full model (multiple linear regression) statistics for alpha diversity in all age
write.table(lm.stats.tbl, file = "../results/tables/alpha_lm_stats.txt", sep = "\t", row.names = F, quote = F)
lm.stats.tbl.2 = lm.stats.tbl %>% dplyr::filter(!is.na(Eta2_partial)) %>% rename(p_value = `Pr(>F)`) %>% mutate(Factors = factor(Factors, levels = c("Age_num", "Gender", "Expo", "Age_num:Expo", "Age_num:Gender", "Expo:Gender"))) %>% arrange(Body_Site, Indice, Year, Factors)
ggplot(lm.stats.tbl.2, aes(x = Factors, y = p_value, color = Body_Site)) +
geom_point() +
geom_path(aes(group = Year, linetype = Year)) +
facet_grid(Body_Site~Indice) +
scale_y_continuous(breaks = c(0, 0.05, 0.1, 0.2, 1.0)) +
coord_cartesian(ylim = c(0, 0.2)) +
theme_jw() + theme(aspect.ratio = 0.6, axis.text.x = element_text(angle = 90, hjust = 1)) +
guides(color = F) +
ggtitle("FDR adjusted p-value for full models of alpha diversity in all age")
ggplot(lm.stats.tbl.2, aes(x = Factors, y = Eta2_partial, color = Body_Site)) +
geom_point() +
geom_path(aes(group = Year, linetype = Year)) +
facet_grid(Body_Site~Indice, scales = "free_y") +
theme_jw() + theme(aspect.ratio = 0.6, axis.text.x = element_text(angle = 90, hjust = 1)) +
guides(color = F) +
ggtitle("Effect size (eta-partial) for full models of alpha diversity in all age")
#p1 + p2 + plot_layout(guides = "collect")Final Models (No interaction terms)
models = list()
model_tbls = data.frame(BodySite = as.character(),
Index = as.character(),
Year = as.character(),
Model = as.character(),
Model.type = as.character())
for (bb in BodySites){
for (aa in Indices){
for (yy in Yrs){
Work_dat <- dat_alpha %>% dplyr::filter(Body_Site==bb, Year==yy, SampleGroup=="Villagers", Age_num>=1)
# Obtain the full linear model
fit.full <- lm(paste0(aa, " ~ Age_num + Expo + Gender"), data = Work_dat)
# Obtain the null model with no predictors, this is necessary for forward selection use stepAIC
fit.null <- lm(paste0(aa, " ~ 1"), data = Work_dat)
# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~ Age_num + Expo + Gender)
lm.model <- lm(regforward.out$terms, data = Work_dat, contrasts = list(Expo = "contr.sum", Gender = "contr.sum"))
lm.model.anova = car::Anova(lm.model, type = 3)
if ("Age_num" %in% rownames(lm.model.anova)){
fit.seg <- segmented(regforward.out, seg.Z = ~ Age_num, psi = 9)
model_comp = anova(fit.seg, lm.model)
if (!is.na(model_comp[2, 6]) & model_comp[2, 6]<=0.05){
fit.final = fit.seg
summary(fit.seg) %>% print()
davies.test(lm.model, seg.Z = ~ Age_num, k = 10) %>% print()
confint.segmented(fit.seg) %>% print()
slope(fit.seg) %>% print()
} else {
fit.final = lm.model
}
} else {
fit.final = lm.model
}
plot(fit.final$fitted.values, fit.final$residuals)
plot((fit.final$model)[, 1], fit.final$fitted.values)
plot(hatvalues(fit.final), fit.final$residuals)
qqnorm(fit.final$residuals)
ff = (fit.final$terms %>% as.character())[3]
models[[bb]][[aa]][[yy]] = fit.final
if ("segmented" %in% class(fit.final)){
model_tbls = rbind(model_tbls, cbind(bb, aa, yy, ff, "segmented"))
} else {
model_tbls = rbind(model_tbls, cbind(bb, aa, yy, ff, "linear"))
}
}
}
}Start: AIC=577.4
faith_pd ~ 1
Df Sum of Sq RSS AIC
+ Expo 1 877.60 5581.9 556.91
+ Age_num 1 251.81 6207.6 573.27
+ Gender 1 126.68 6332.8 576.35
<none> 6459.5 577.40
Step: AIC=556.91
faith_pd ~ Expo
Df Sum of Sq RSS AIC
+ Gender 1 177.43 5404.4 553.94
+ Age_num 1 134.92 5446.9 555.14
<none> 5581.9 556.91
Step: AIC=553.94
faith_pd ~ Expo + Gender
Df Sum of Sq RSS AIC
+ Age_num 1 164.54 5239.9 551.17
<none> 5404.4 553.94
Step: AIC=551.17
faith_pd ~ Expo + Gender + Age_num
***Regression Model with Segmented Relationship(s)***
Call:
segmented.lm(obj = regforward.out, seg.Z = ~Age_num, psi = 9)
Estimated Break-Point(s):
Est. St.Err
psi1.Age_num 7.918 1.69
Meaningful coefficients of the linear terms:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 20.2457 2.0376 9.936 < 2e-16 ***
ExpoMedium -4.3660 0.9399 -4.645 7.45e-06 ***
GenderM -2.2691 0.9107 -2.492 0.01382 *
Age_num 1.1480 0.4147 2.768 0.00636 **
U1.Age_num -1.1658 0.4157 -2.804 NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 5.578 on 148 degrees of freedom
Multiple R-Squared: 0.287, Adjusted R-squared: 0.2629
Convergence attained in 3 iter. (rel. change 1.9747e-16)
Davies' test for a change in the slope
data: formula = faith_pd ~ Expo + Gender + Age_num , method = lm
model = gaussian , link = identity
segmented variable = Age_num
'best' at = 7.5556, n.points = 9, p-value = 0.0001387
alternative hypothesis: two.sided
Est. CI(95%).low CI(95%).up
psi1.Age_num 7.91839 4.57902 11.2578
$Age_num
Est. St.Err. t value CI(95%).l CI(95%).u
slope1 1.148000 0.414700 2.76820 0.328470 1.967500
slope2 -0.017832 0.035963 -0.49584 -0.088899 0.053236
Start: AIC=610.01
faith_pd ~ 1
Df Sum of Sq RSS AIC
+ Expo 1 466.10 5026.9 596.32
+ Age_num 1 74.17 5418.9 609.60
<none> 5493.0 610.01
+ Gender 1 34.72 5458.3 610.89
Step: AIC=596.32
faith_pd ~ Expo
Df Sum of Sq RSS AIC
<none> 5026.9 596.32
+ Gender 1 55.767 4971.2 596.34
+ Age_num 1 44.091 4982.8 596.76
Start: AIC=-798.63
pielou_e ~ 1
Df Sum of Sq RSS AIC
<none> 0.85047 -798.63
+ Expo 1 0.0072315 0.84324 -797.95
+ Gender 1 0.0067885 0.84368 -797.87
+ Age_num 1 0.0013146 0.84915 -796.87
Start: AIC=-928.53
pielou_e ~ 1
Df Sum of Sq RSS AIC
+ Expo 1 0.048001 0.87411 -935.99
<none> 0.92211 -928.53
+ Age_num 1 0.000622 0.92149 -926.65
+ Gender 1 0.000009 0.92210 -926.53
Step: AIC=-935.99
pielou_e ~ Expo
Df Sum of Sq RSS AIC
<none> 0.87411 -935.99
+ Gender 1 0.00035701 0.87375 -934.07
+ Age_num 1 0.00002161 0.87409 -934.00
Start: AIC=-52.93
shannon ~ 1
Df Sum of Sq RSS AIC
+ Expo 1 5.8680 101.93 -59.555
+ Gender 1 1.9180 105.88 -53.699
<none> 107.80 -52.935
+ Age_num 1 0.7776 107.02 -52.050
Step: AIC=-59.55
shannon ~ Expo
Df Sum of Sq RSS AIC
+ Gender 1 2.41770 99.51 -61.251
<none> 101.93 -59.555
+ Age_num 1 0.28031 101.65 -57.979
Step: AIC=-61.25
shannon ~ Expo + Gender
Df Sum of Sq RSS AIC
<none> 99.510 -61.251
+ Age_num 1 0.44589 99.064 -59.943
Start: AIC=-99.94
shannon ~ 1
Df Sum of Sq RSS AIC
+ Expo 1 7.3360 92.171 -111.493
<none> 99.507 -99.938
+ Age_num 1 0.3042 99.203 -98.480
+ Gender 1 0.0459 99.461 -98.019
Step: AIC=-111.49
shannon ~ Expo
Df Sum of Sq RSS AIC
<none> 92.171 -111.49
+ Gender 1 0.168447 92.003 -109.82
+ Age_num 1 0.091109 92.080 -109.67
Start: AIC=1366.98
observed_otus ~ 1
Df Sum of Sq RSS AIC
+ Expo 1 109046 979633 1352.7
+ Age_num 1 41607 1047072 1363.0
+ Gender 1 21417 1067263 1365.9
<none> 1088679 1367.0
Step: AIC=1352.73
observed_otus ~ Expo
Df Sum of Sq RSS AIC
+ Gender 1 28701 950932 1350.2
+ Age_num 1 24608 955025 1350.8
<none> 979633 1352.7
Step: AIC=1350.15
observed_otus ~ Expo + Gender
Df Sum of Sq RSS AIC
+ Age_num 1 29694 921238 1347.3
<none> 950932 1350.2
Step: AIC=1347.26
observed_otus ~ Expo + Gender + Age_num
***Regression Model with Segmented Relationship(s)***
Call:
segmented.lm(obj = regforward.out, seg.Z = ~Age_num, psi = 9)
Estimated Break-Point(s):
Est. St.Err
psi1.Age_num 8.333 1.747
Meaningful coefficients of the linear terms:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 170.512 24.597 6.932 1.19e-10 ***
ExpoMedium -48.451 12.635 -3.835 0.000186 ***
GenderM -29.088 12.246 -2.375 0.018818 *
Age_num 12.481 4.284 2.913 0.004131 **
U1.Age_num -12.605 4.307 -2.927 NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 75.27 on 148 degrees of freedom
Multiple R-Squared: 0.2298, Adjusted R-squared: 0.2038
Convergence attained in 4 iter. (rel. change 0)
Davies' test for a change in the slope
data: formula = observed_otus ~ Expo + Gender + Age_num , method = lm
model = gaussian , link = identity
segmented variable = Age_num
'best' at = 7.5556, n.points = 9, p-value = 0.001613
alternative hypothesis: two.sided
Est. CI(95%).low CI(95%).up
psi1.Age_num 8.3333 4.88077 11.7858
$Age_num
Est. St.Err. t value CI(95%).l CI(95%).u
slope1 12.48100 4.28410 2.91330 4.0150 20.94700
slope2 -0.12465 0.52149 -0.23902 -1.1552 0.90588
Start: AIC=1452.83
observed_otus ~ 1
Df Sum of Sq RSS AIC
+ Expo 1 34581 607797 1445.0
<none> 642379 1452.8
+ Age_num 1 3124 639254 1454.0
+ Gender 1 1626 640753 1454.4
Step: AIC=1445.04
observed_otus ~ Expo
Df Sum of Sq RSS AIC
<none> 607797 1445.0
+ Gender 1 2901.4 604896 1446.2
+ Age_num 1 1507.6 606290 1446.6
Start: AIC=503.39
faith_pd ~ 1
Df Sum of Sq RSS AIC
+ Expo 1 327.20 3160.9 489.23
+ Age_num 1 326.20 3161.9 489.29
<none> 3488.1 503.39
+ Gender 1 3.18 3484.9 505.24
Step: AIC=489.23
faith_pd ~ Expo
Df Sum of Sq RSS AIC
+ Age_num 1 218.027 2942.9 479.51
<none> 3160.9 489.23
+ Gender 1 0.781 3160.1 491.19
Step: AIC=479.51
faith_pd ~ Expo + Age_num
Df Sum of Sq RSS AIC
<none> 2942.9 479.51
+ Gender 1 0.062452 2942.8 481.51
Start: AIC=448.86
faith_pd ~ 1
Df Sum of Sq RSS AIC
+ Age_num 1 578.15 2494.4 420.21
+ Expo 1 224.27 2848.3 439.71
<none> 3072.6 448.86
+ Gender 1 31.84 3040.7 449.32
Step: AIC=420.21
faith_pd ~ Age_num
Df Sum of Sq RSS AIC
+ Expo 1 172.536 2321.9 411.68
<none> 2494.4 420.21
+ Gender 1 3.114 2491.3 422.03
Step: AIC=411.68
faith_pd ~ Age_num + Expo
Df Sum of Sq RSS AIC
<none> 2321.9 411.68
+ Gender 1 0.87557 2321.0 413.62
Start: AIC=-736.73
pielou_e ~ 1
Df Sum of Sq RSS AIC
+ Expo 1 0.33382 1.4799 -768.09
<none> 1.8138 -736.73
+ Gender 1 0.01466 1.7991 -736.06
+ Age_num 1 0.00116 1.8126 -734.84
Step: AIC=-768.09
pielou_e ~ Expo
Df Sum of Sq RSS AIC
+ Gender 1 0.0225098 1.4574 -768.60
<none> 1.4799 -768.09
+ Age_num 1 0.0067515 1.4732 -766.84
Step: AIC=-768.6
pielou_e ~ Expo + Gender
Df Sum of Sq RSS AIC
<none> 1.4574 -768.60
+ Age_num 1 0.0050221 1.4524 -767.17
Start: AIC=-674.18
pielou_e ~ 1
Df Sum of Sq RSS AIC
<none> 1.4779 -674.18
+ Expo 1 0.0086101 1.4692 -673.04
+ Gender 1 0.0045496 1.4733 -672.63
+ Age_num 1 0.0041333 1.4737 -672.59
Start: AIC=-16.6
shannon ~ 1
Df Sum of Sq RSS AIC
+ Expo 1 30.7103 115.70 -53.211
<none> 146.41 -16.604
+ Age_num 1 1.1004 145.31 -15.841
+ Gender 1 0.3191 146.09 -14.962
Step: AIC=-53.21
shannon ~ Expo
Df Sum of Sq RSS AIC
<none> 115.70 -53.211
+ Gender 1 0.70883 114.99 -52.219
+ Age_num 1 0.00265 115.70 -51.215
Start: AIC=-10.31
shannon ~ 1
Df Sum of Sq RSS AIC
+ Age_num 1 3.6509 131.54 -12.3362
<none> 135.19 -10.3118
+ Gender 1 0.4367 134.75 -8.7875
+ Expo 1 0.2287 134.96 -8.5607
Step: AIC=-12.34
shannon ~ Age_num
Df Sum of Sq RSS AIC
<none> 131.54 -12.336
+ Gender 1 0.12639 131.41 -10.478
+ Expo 1 0.10877 131.43 -10.458
Start: AIC=1202.36
observed_otus ~ 1
Df Sum of Sq RSS AIC
+ Expo 1 25991.7 221496 1186.2
+ Age_num 1 9449.1 238038 1198.0
<none> 247488 1202.4
+ Gender 1 191.2 247296 1204.2
Step: AIC=1186.16
observed_otus ~ Expo
Df Sum of Sq RSS AIC
+ Age_num 1 4428.0 217068 1184.8
<none> 221496 1186.2
+ Gender 1 33.7 221462 1188.1
Step: AIC=1184.85
observed_otus ~ Expo + Age_num
Df Sum of Sq RSS AIC
<none> 217068 1184.8
+ Gender 1 0.49699 217067 1186.8
Start: AIC=1099.97
observed_otus ~ 1
Df Sum of Sq RSS AIC
+ Age_num 1 28638.6 229086 1084.7
+ Expo 1 13366.0 244359 1094.2
<none> 257725 1100.0
+ Gender 1 416.3 257309 1101.7
Step: AIC=1084.66
observed_otus ~ Age_num
Df Sum of Sq RSS AIC
+ Expo 1 10545.4 218541 1079.7
<none> 229086 1084.7
+ Gender 1 51.1 229035 1086.6
Step: AIC=1079.73
observed_otus ~ Age_num + Expo
Df Sum of Sq RSS AIC
<none> 218541 1079.7
+ Gender 1 186.96 218354 1081.6
Start: AIC=717.39
faith_pd ~ 1
Df Sum of Sq RSS AIC
+ Expo 1 1204.32 24178 712.73
+ Age_num 1 784.55 24598 715.09
<none> 25383 717.39
+ Gender 1 1.37 25381 719.39
Step: AIC=712.73
faith_pd ~ Expo
Df Sum of Sq RSS AIC
+ Age_num 1 1193.82 22985 707.80
<none> 24178 712.73
+ Gender 1 0.87 24178 714.73
Step: AIC=707.8
faith_pd ~ Expo + Age_num
Df Sum of Sq RSS AIC
<none> 22985 707.80
+ Gender 1 13.127 22972 709.72
***Regression Model with Segmented Relationship(s)***
Call:
segmented.lm(obj = regforward.out, seg.Z = ~Age_num, psi = 9)
Estimated Break-Point(s):
Est. St.Err
psi1.Age_num 10.503 1.915
Meaningful coefficients of the linear terms:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.2593 3.5611 0.354 0.724183
ExpoMedium 7.4982 2.1762 3.446 0.000764 ***
Age_num 2.0584 0.5730 3.592 0.000461 ***
U1.Age_num -2.1586 0.5797 -3.724 NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 12.12 on 132 degrees of freedom
Multiple R-Squared: 0.2359, Adjusted R-squared: 0.2127
Convergence attained in 3 iter. (rel. change 0)
Davies' test for a change in the slope
data: formula = faith_pd ~ Expo + Age_num , method = lm
model = gaussian , link = identity
segmented variable = Age_num
'best' at = 14.111, n.points = 9, p-value = 7.96e-05
alternative hypothesis: two.sided
Est. CI(95%).low CI(95%).up
psi1.Age_num 10.5028 6.71477 14.2907
$Age_num
Est. St.Err. t value CI(95%).l CI(95%).u
slope1 2.05840 0.573020 3.5922 0.92492 3.191900
slope2 -0.10021 0.096118 -1.0426 -0.29034 0.089921
Start: AIC=661.41
faith_pd ~ 1
Df Sum of Sq RSS AIC
+ Age_num 1 797.06 15612 656.54
+ Expo 1 287.61 16122 660.97
<none> 16409 661.41
+ Gender 1 14.96 16394 663.29
Step: AIC=656.54
faith_pd ~ Age_num
Df Sum of Sq RSS AIC
+ Expo 1 515.66 15097 653.91
<none> 15612 656.54
+ Gender 1 82.15 15530 657.81
Step: AIC=653.91
faith_pd ~ Age_num + Expo
Df Sum of Sq RSS AIC
<none> 15097 653.91
+ Gender 1 64.915 15032 655.31
***Regression Model with Segmented Relationship(s)***
Call:
segmented.lm(obj = regforward.out, seg.Z = ~Age_num, psi = 9)
Estimated Break-Point(s):
Est. St.Err
psi1.Age_num 13 2.365
Meaningful coefficients of the linear terms:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.3916 2.8886 1.867 0.064172 .
Age_num 1.3496 0.3750 3.599 0.000449 ***
ExpoMedium 4.6226 1.8022 2.565 0.011427 *
U1.Age_num -1.4254 0.3823 -3.729 NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 9.854 on 133 degrees of freedom
Multiple R-Squared: 0.213, Adjusted R-squared: 0.1894
Convergence attained in 5 iter. (rel. change 3.6307e-07)
Davies' test for a change in the slope
data: formula = faith_pd ~ Age_num + Expo , method = lm
model = gaussian , link = identity
segmented variable = Age_num
'best' at = 15.444, n.points = 9, p-value = 0.0002563
alternative hypothesis: two.sided
Est. CI(95%).low CI(95%).up
psi1.Age_num 12.9999 8.32219 17.6777
$Age_num
Est. St.Err. t value CI(95%).l CI(95%).u
slope1 1.349600 0.374950 3.59950 0.60799 2.091300
slope2 -0.075799 0.082967 -0.91361 -0.23990 0.088305
Start: AIC=-563.47
pielou_e ~ 1
Df Sum of Sq RSS AIC
+ Expo 1 0.040152 2.1687 -563.98
<none> 2.2088 -563.47
+ Age_num 1 0.017469 2.1914 -562.56
+ Gender 1 0.003786 2.2050 -561.70
Step: AIC=-563.98
pielou_e ~ Expo
Df Sum of Sq RSS AIC
<none> 2.1687 -563.98
+ Age_num 1 0.0097853 2.1589 -562.60
+ Gender 1 0.0036188 2.1651 -562.21
Start: AIC=-586.78
pielou_e ~ 1
Df Sum of Sq RSS AIC
+ Expo 1 0.043216 1.8930 -587.89
<none> 1.9363 -586.78
+ Gender 1 0.008933 1.9273 -585.41
+ Age_num 1 0.006969 1.9293 -585.27
Step: AIC=-587.89
pielou_e ~ Expo
Df Sum of Sq RSS AIC
<none> 1.8930 -587.89
+ Gender 1 0.0122879 1.8808 -586.79
+ Age_num 1 0.0020227 1.8910 -586.04
Start: AIC=67.68
shannon ~ 1
Df Sum of Sq RSS AIC
+ Age_num 1 11.5234 209.75 62.353
<none> 221.27 67.680
+ Gender 1 0.2747 221.00 69.510
+ Expo 1 0.2379 221.03 69.533
Step: AIC=62.35
shannon ~ Age_num
Df Sum of Sq RSS AIC
<none> 209.75 62.353
+ Gender 1 0.93295 208.82 63.742
+ Expo 1 0.01050 209.74 64.346
***Regression Model with Segmented Relationship(s)***
Call:
segmented.lm(obj = regforward.out, seg.Z = ~Age_num, psi = 9)
Estimated Break-Point(s):
Est. St.Err
psi1.Age_num 12.022 2.763
Meaningful coefficients of the linear terms:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.94666 0.31238 6.232 5.66e-09 ***
Age_num 0.14656 0.05043 2.906 0.00429 **
U1.Age_num -0.15464 0.05135 -3.012 NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.182 on 133 degrees of freedom
Multiple R-Squared: 0.1598, Adjusted R-squared: 0.1409
Convergence attained in 3 iter. (rel. change 3.0031e-07)
Davies' test for a change in the slope
data: formula = shannon ~ Age_num , method = lm
model = gaussian , link = identity
segmented variable = Age_num
'best' at = 14.111, n.points = 9, p-value = 0.0006814
alternative hypothesis: two.sided
Est. CI(95%).low CI(95%).up
psi1.Age_num 12.0223 6.55758 17.487
$Age_num
Est. St.Err. t value CI(95%).l CI(95%).u
slope1 0.1465600 0.0504330 2.90610 0.046809 0.246320
slope2 -0.0080792 0.0096185 -0.83996 -0.027104 0.010946
Start: AIC=40.36
shannon ~ 1
Df Sum of Sq RSS AIC
+ Age_num 1 7.5644 174.65 36.506
<none> 182.22 40.357
+ Expo 1 0.4060 181.81 42.049
+ Gender 1 0.1408 182.08 42.250
Step: AIC=36.51
shannon ~ Age_num
Df Sum of Sq RSS AIC
<none> 174.65 36.506
+ Gender 1 0.77672 173.88 37.891
+ Expo 1 0.01420 174.64 38.494
***Regression Model with Segmented Relationship(s)***
Call:
segmented.lm(obj = regforward.out, seg.Z = ~Age_num, psi = 9)
Estimated Break-Point(s):
Est. St.Err
psi1.Age_num 17.227 3.354
Meaningful coefficients of the linear terms:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.17025 0.24915 8.710 9.96e-15 ***
Age_num 0.09580 0.02856 3.355 0.00103 **
U1.Age_num -0.11278 0.03039 -3.711 NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.069 on 134 degrees of freedom
Multiple R-Squared: 0.1591, Adjusted R-squared: 0.1403
Convergence attained in 4 iter. (rel. change 0)
Davies' test for a change in the slope
data: formula = shannon ~ Age_num , method = lm
model = gaussian , link = identity
segmented variable = Age_num
'best' at = 15.444, n.points = 9, p-value = 0.0002751
alternative hypothesis: two.sided
Est. CI(95%).low CI(95%).up
psi1.Age_num 17.2273 10.5927 23.8619
$Age_num
Est. St.Err. t value CI(95%).l CI(95%).u
slope1 0.095802 0.028556 3.3549 0.039324 0.1522800
slope2 -0.016981 0.010403 -1.6322 -0.037557 0.0035952
Start: AIC=1264.61
observed_otus ~ 1
Df Sum of Sq RSS AIC
+ Expo 1 45099 1332819 1262.0
<none> 1377918 1264.6
+ Age_num 1 19166 1358752 1264.7
+ Gender 1 758 1377160 1266.5
Step: AIC=1262.05
observed_otus ~ Expo
Df Sum of Sq RSS AIC
+ Age_num 1 31672 1301147 1260.8
<none> 1332819 1262.0
+ Gender 1 840 1331979 1264.0
Step: AIC=1260.75
observed_otus ~ Expo + Age_num
Df Sum of Sq RSS AIC
<none> 1301147 1260.8
+ Gender 1 2782.7 1298365 1262.5
***Regression Model with Segmented Relationship(s)***
Call:
segmented.lm(obj = regforward.out, seg.Z = ~Age_num, psi = 9)
Estimated Break-Point(s):
Est. St.Err
psi1.Age_num 10 1.86
Meaningful coefficients of the linear terms:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -14.592 27.044 -0.540 0.590405
ExpoMedium 45.538 16.526 2.755 0.006690 **
Age_num 15.130 4.352 3.477 0.000687 ***
U1.Age_num -16.118 4.402 -3.661 NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 92.05 on 132 degrees of freedom
Multiple R-Squared: 0.1882, Adjusted R-squared: 0.1636
Convergence attained in 6 iter. (rel. change 1.8623e-07)
Davies' test for a change in the slope
data: formula = observed_otus ~ Expo + Age_num , method = lm
model = gaussian , link = identity
segmented variable = Age_num
'best' at = 7.5556, n.points = 9, p-value = 0.0002517
alternative hypothesis: two.sided
Est. CI(95%).low CI(95%).up
psi1.Age_num 10 6.3204 13.6796
$Age_num
Est. St.Err. t value CI(95%).l CI(95%).u
slope1 15.130 4.35170 3.4769 6.5225 23.7380
slope2 -0.988 0.72994 -1.3535 -2.4319 0.4559
Start: AIC=1188.48
observed_otus ~ 1
Df Sum of Sq RSS AIC
+ Age_num 1 24477.4 723352 1185.9
<none> 747829 1188.5
+ Expo 1 4861.7 742967 1189.6
+ Gender 1 1190.4 746639 1190.3
Step: AIC=1185.89
observed_otus ~ Age_num
Df Sum of Sq RSS AIC
<none> 723352 1185.9
+ Expo 1 10229 713122 1185.9
+ Gender 1 4033 719319 1187.1
***Regression Model with Segmented Relationship(s)***
Call:
segmented.lm(obj = regforward.out, seg.Z = ~Age_num, psi = 9)
Estimated Break-Point(s):
Est. St.Err
psi1.Age_num 13 2.513
Meaningful coefficients of the linear terms:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 28.101 17.880 1.572 0.118390
Age_num 8.146 2.338 3.485 0.000666 ***
U1.Age_num -8.876 2.416 -3.675 NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 69.02 on 134 degrees of freedom
Multiple R-Squared: 0.1463, Adjusted R-squared: 0.1272
Convergence attained in 6 iter. (rel. change 7.0875e-06)
Davies' test for a change in the slope
data: formula = observed_otus ~ Age_num , method = lm
model = gaussian , link = identity
segmented variable = Age_num
'best' at = 15.444, n.points = 9, p-value = 0.0009512
alternative hypothesis: two.sided
Est. CI(95%).low CI(95%).up
psi1.Age_num 13.0001 8.03041 17.9698
$Age_num
Est. St.Err. t value CI(95%).l CI(95%).u
slope1 8.14620 2.33770 3.4848 3.5227 12.77000
slope2 -0.73013 0.60844 -1.2000 -1.9335 0.47326
Start: AIC=898.45
faith_pd ~ 1
Df Sum of Sq RSS AIC
+ Expo 1 9453.0 29340 854.64
+ Gender 1 750.3 38043 897.24
<none> 38793 898.45
+ Age_num 1 301.5 38492 899.17
Step: AIC=854.64
faith_pd ~ Expo
Df Sum of Sq RSS AIC
+ Gender 1 978.22 28362 851.08
<none> 29340 854.64
+ Age_num 1 0.12 29340 856.64
Step: AIC=851.08
faith_pd ~ Expo + Gender
Df Sum of Sq RSS AIC
<none> 28362 851.08
+ Age_num 1 1.1997 28361 853.07
Start: AIC=457.7
faith_pd ~ 1
Df Sum of Sq RSS AIC
+ Expo 1 1468.41 18646 453.41
<none> 20115 457.70
+ Gender 1 61.22 20054 459.45
+ Age_num 1 51.33 20064 459.49
Step: AIC=453.41
faith_pd ~ Expo
Df Sum of Sq RSS AIC
<none> 18646 453.41
+ Age_num 1 107.687 18539 454.93
+ Gender 1 52.232 18594 455.18
Start: AIC=-905.64
pielou_e ~ 1
Df Sum of Sq RSS AIC
+ Gender 1 0.0292817 0.61829 -911.23
<none> 0.64757 -905.64
+ Expo 1 0.0052636 0.64231 -904.98
+ Age_num 1 0.0004897 0.64708 -903.76
Step: AIC=-911.23
pielou_e ~ Gender
Df Sum of Sq RSS AIC
<none> 0.61829 -911.23
+ Expo 1 0.0063057 0.61199 -910.91
+ Age_num 1 0.0009688 0.61732 -909.49
Start: AIC=-453.82
pielou_e ~ 1
Df Sum of Sq RSS AIC
+ Expo 1 0.040233 0.30175 -462.21
+ Age_num 1 0.011223 0.33076 -454.59
<none> 0.34199 -453.82
+ Gender 1 0.001064 0.34092 -452.08
Step: AIC=-462.21
pielou_e ~ Expo
Df Sum of Sq RSS AIC
+ Age_num 1 0.0151271 0.28663 -464.48
<none> 0.30175 -462.21
+ Gender 1 0.0012784 0.30047 -460.56
Step: AIC=-464.48
pielou_e ~ Expo + Age_num
Df Sum of Sq RSS AIC
<none> 0.28663 -464.48
+ Gender 1 0.0016638 0.28496 -462.96
Start: AIC=-67.4
shannon ~ 1
Df Sum of Sq RSS AIC
+ Expo 1 6.6969 100.72 -75.954
+ Gender 1 4.8961 102.52 -73.048
<none> 107.42 -67.397
+ Age_num 1 0.2545 107.16 -65.786
Step: AIC=-75.95
shannon ~ Expo
Df Sum of Sq RSS AIC
+ Gender 1 5.3699 95.35 -82.940
<none> 100.72 -75.954
+ Age_num 1 0.0012 100.72 -73.956
Step: AIC=-82.94
shannon ~ Expo + Gender
Df Sum of Sq RSS AIC
<none> 95.35 -82.940
+ Age_num 1 0.019726 95.33 -80.974
Start: AIC=-33.65
shannon ~ 1
Df Sum of Sq RSS AIC
+ Expo 1 6.8840 47.136 -42.962
<none> 54.020 -33.647
+ Age_num 1 1.0369 52.983 -33.256
+ Gender 1 0.0162 54.004 -31.672
Step: AIC=-42.96
shannon ~ Expo
Df Sum of Sq RSS AIC
+ Age_num 1 1.53787 45.598 -43.715
<none> 47.136 -42.962
+ Gender 1 0.02836 47.108 -41.012
Step: AIC=-43.71
shannon ~ Expo + Age_num
Df Sum of Sq RSS AIC
<none> 45.598 -43.715
+ Gender 1 0.04796 45.550 -41.802
Start: AIC=1575.04
observed_otus ~ 1
Df Sum of Sq RSS AIC
+ Expo 1 415511 1985849 1545.9
+ Gender 1 58732 2342629 1573.0
<none> 2401361 1575.0
+ Age_num 1 13739 2387621 1576.1
Step: AIC=1545.88
observed_otus ~ Expo
Df Sum of Sq RSS AIC
+ Gender 1 71908 1913942 1541.8
<none> 1985849 1545.9
+ Age_num 1 0 1985849 1547.9
Step: AIC=1541.83
observed_otus ~ Expo + Gender
Df Sum of Sq RSS AIC
<none> 1913942 1541.8
+ Age_num 1 148.46 1913793 1543.8
Start: AIC=822.68
observed_otus ~ 1
Df Sum of Sq RSS AIC
+ Expo 1 62837 1571116 821.42
<none> 1633952 822.68
+ Gender 1 4528 1629424 824.45
+ Age_num 1 3480 1630472 824.50
Step: AIC=821.42
observed_otus ~ Expo
Df Sum of Sq RSS AIC
<none> 1571116 821.42
+ Age_num 1 6408.9 1564707 823.08
+ Gender 1 4018.0 1567098 823.21
Start: AIC=912.08
faith_pd ~ 1
Df Sum of Sq RSS AIC
+ Expo 1 17233.1 24922 827.88
+ Age_num 1 1860.3 40295 906.67
+ Gender 1 1403.9 40751 908.52
<none> 42155 912.08
Step: AIC=827.88
faith_pd ~ Expo
Df Sum of Sq RSS AIC
+ Gender 1 1705.46 23217 818.25
+ Age_num 1 532.68 24389 826.33
<none> 24922 827.88
Step: AIC=818.25
faith_pd ~ Expo + Gender
Df Sum of Sq RSS AIC
+ Age_num 1 588.06 22629 816.04
<none> 23217 818.25
Step: AIC=816.04
faith_pd ~ Expo + Gender + Age_num
Start: AIC=665.58
faith_pd ~ 1
Df Sum of Sq RSS AIC
+ Expo 1 7975.2 26015 636.30
+ Age_num 1 789.1 33201 664.84
+ Gender 1 641.0 33349 665.36
<none> 33990 665.58
Step: AIC=636.3
faith_pd ~ Expo
Df Sum of Sq RSS AIC
+ Gender 1 799.88 25215 634.64
+ Age_num 1 536.00 25479 635.86
<none> 26015 636.30
Step: AIC=634.64
faith_pd ~ Expo + Gender
Df Sum of Sq RSS AIC
+ Age_num 1 545.27 24670 634.09
<none> 25215 634.64
Step: AIC=634.09
faith_pd ~ Expo + Gender + Age_num
Start: AIC=-916.45
pielou_e ~ 1
Df Sum of Sq RSS AIC
+ Expo 1 0.0310501 0.57520 -923.08
+ Age_num 1 0.0218545 0.58440 -920.47
<none> 0.60625 -916.45
+ Gender 1 0.0044785 0.60178 -915.67
Step: AIC=-923.08
pielou_e ~ Expo
Df Sum of Sq RSS AIC
+ Age_num 1 0.0148881 0.56032 -925.38
<none> 0.57520 -923.08
+ Gender 1 0.0051937 0.57001 -922.56
Step: AIC=-925.38
pielou_e ~ Expo + Age_num
Df Sum of Sq RSS AIC
<none> 0.56032 -925.38
+ Gender 1 0.0057042 0.55461 -925.05
Start: AIC=-538.22
pielou_e ~ 1
Df Sum of Sq RSS AIC
+ Age_num 1 0.066373 1.0896 -543.14
+ Expo 1 0.020103 1.1358 -538.27
<none> 1.1559 -538.22
+ Gender 1 0.007772 1.1482 -537.01
Step: AIC=-543.14
pielou_e ~ Age_num
Df Sum of Sq RSS AIC
<none> 1.0896 -543.14
+ Expo 1 0.0162888 1.0733 -542.90
+ Gender 1 0.0081806 1.0814 -542.02
***Regression Model with Segmented Relationship(s)***
Call:
segmented.lm(obj = regforward.out, seg.Z = ~Age_num, psi = 9)
Estimated Break-Point(s):
Est. St.Err
psi1.Age_num 44.999 6.518
Meaningful coefficients of the linear terms:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.8078937 0.0153066 52.781 < 2e-16 ***
Age_num -0.0027632 0.0008392 -3.293 0.00132 **
U1.Age_num 0.0085514 0.0041917 2.040 NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.0955 on 113 degrees of freedom
Multiple R-Squared: 0.1085, Adjusted R-squared: 0.08479
Convergence attained in 1 iter. (rel. change 3.4462e-07)
Davies' test for a change in the slope
data: formula = pielou_e ~ Age_num , method = lm
model = gaussian , link = identity
segmented variable = Age_num
'best' at = 46.889, n.points = 9, p-value = 0.09208
alternative hypothesis: two.sided
Est. CI(95%).low CI(95%).up
psi1.Age_num 44.9993 32.0859 57.9127
$Age_num
Est. St.Err. t value CI(95%).l CI(95%).u
slope1 -0.0027632 0.00083918 -3.2927 -0.0044258 -0.0011006
slope2 0.0057882 0.00410680 1.4094 -0.0023482 0.0139250
Start: AIC=-76.15
shannon ~ 1
Df Sum of Sq RSS AIC
+ Expo 1 19.5685 82.264 -109.149
+ Age_num 1 4.6701 97.162 -81.851
+ Gender 1 2.0289 99.804 -77.452
<none> 101.833 -76.152
Step: AIC=-109.15
shannon ~ Expo
Df Sum of Sq RSS AIC
+ Gender 1 2.4132 79.851 -112.03
+ Age_num 1 2.2321 80.032 -111.66
<none> 82.264 -109.15
Step: AIC=-112.03
shannon ~ Expo + Gender
Df Sum of Sq RSS AIC
+ Age_num 1 2.3663 77.484 -114.97
<none> 79.851 -112.03
Step: AIC=-114.97
shannon ~ Expo + Gender + Age_num
Start: AIC=33.32
shannon ~ 1
Df Sum of Sq RSS AIC
+ Expo 1 13.5753 139.34 24.446
+ Age_num 1 8.3542 144.56 28.750
<none> 152.92 33.323
+ Gender 1 2.1475 150.77 33.668
Step: AIC=24.45
shannon ~ Expo
Df Sum of Sq RSS AIC
+ Age_num 1 7.2321 132.11 20.210
+ Gender 1 2.5217 136.82 24.309
<none> 139.34 24.446
Step: AIC=20.21
shannon ~ Expo + Age_num
Df Sum of Sq RSS AIC
+ Gender 1 2.5822 129.53 19.901
<none> 132.11 20.210
Step: AIC=19.9
shannon ~ Expo + Age_num + Gender
Start: AIC=1592.59
observed_otus ~ 1
Df Sum of Sq RSS AIC
+ Expo 1 938567 1734140 1523.7
+ Age_num 1 106387 2566320 1587.9
+ Gender 1 62096 2610611 1590.7
<none> 2672707 1592.6
Step: AIC=1523.65
observed_otus ~ Expo
Df Sum of Sq RSS AIC
+ Gender 1 76969 1657171 1518.2
+ Age_num 1 31788 1702352 1522.6
<none> 1734140 1523.7
Step: AIC=1518.2
observed_otus ~ Expo + Gender
Df Sum of Sq RSS AIC
+ Age_num 1 34657 1622514 1516.7
<none> 1657171 1518.2
Step: AIC=1516.74
observed_otus ~ Expo + Gender + Age_num
Start: AIC=1182.19
observed_otus ~ 1
Df Sum of Sq RSS AIC
+ Expo 1 534368 2277124 1159.5
+ Age_num 1 71038 2740454 1181.2
+ Gender 1 57098 2754394 1181.8
<none> 2811492 1182.2
Step: AIC=1159.52
observed_otus ~ Expo
Df Sum of Sq RSS AIC
+ Gender 1 69292 2207832 1157.9
+ Age_num 1 51145 2225979 1158.9
<none> 2277124 1159.5
Step: AIC=1157.91
observed_otus ~ Expo + Gender
Df Sum of Sq RSS AIC
+ Age_num 1 51988 2155844 1157.1
<none> 2207832 1157.9
Step: AIC=1157.12
observed_otus ~ Expo + Gender + Age_num
colnames(model_tbls) = c("Body_Site", "Indice", "Year", "Model", "Model.type")
model_tbls$Model <- gsub(" + psi1.Age_num", "", model_tbls$Model, fixed = T) %>% paste0(model_tbls$Indice, " ~ ", .)Graphs
Indices <- colnames(alpha[c(2:5)])
Yrs <- c("2015", "2016")
BodySites <- c("Feces", "Nose", "Mouth", "Right_Arm", "Right_Hand")
dat_new = complete(cbind(Age_num = seq(1, 60, 1),
Expo = c("Low", "Medium"),
Gender = c("F", "M")) %>%
as.data.frame(), Age_num, Expo, Gender) %>%
mutate(Age_num = as.numeric(Age_num))
dat_pred = list()
for (bb in BodySites){
for (aa in Indices){
for(yy in Yrs){
dat_pred[[bb]][[aa]][[yy]] = cbind(dat_new, predict(models[[bb]][[aa]][[yy]], newdata = dat_new, interval = "confidence"))
}
}
}
for (yy in Yrs){
for (aa in Indices){
#aa = "faith_pd"
p_list = list()
for (bb in BodySites){
#bb = "Feces"
dat_ = dat_alpha %>% dplyr::filter(Body_Site==bb, Year==yy, SampleGroup=="Villagers", Age_num>=1) %>% mutate(Expo = as.character(Expo))
dat_pred_ = dat_pred[[bb]][[aa]][[yy]]
mod = model_tbls %>% dplyr::filter(Body_Site==bb, Year==yy, Indice==aa) %>% pull(Model)
p.expo = (models[[bb]][[aa]][[yy]] %>% car::Anova(type = "3") %>% as.data.frame())["Expo", "Pr(>F)"]
p.expo = ifelse(p.expo<0.001, "P.Expo<0.001", paste0("P.Expo=", round(p.expo, 3)))
p = ggplot(dat_, aes_string(x = "Age_num", y = aa, color = "Expo", shape = "Gender")) +
geom_point(size = 1.5, alpha = 0.6) +
geom_line(data = dat_pred_, aes(y = fit, linetype = Gender), size = 1) +
geom_ribbon(data = dat_pred_, aes(y = fit, ymin = lwr, ymax = upr, fill = Expo, alpha = Gender), color = NA, show.legend = F) +
annotate("label", x = Inf, y = Inf, label = paste(stringi::stri_wrap(c(mod, p.expo),width = 50), collapse = "\n"), hjust = 1, vjust = 1, size = 9/.pt) +
scale_color_manual(values = Expo_color, limits = force) +
scale_fill_manual(values = Expo_color, limits = force) +
scale_shape_manual(values = c(16, 1)) +
scale_linetype_manual(values = c(1, 2)) +
scale_alpha_manual(values = c(0.4, 0.2)) +
theme_jw(base_size = 14) + theme(aspect.ratio = 0.7, panel.background = element_rect(color = "black"), plot.title = element_text(hjust = 0, vjust = 1), plot.title.position = "panel", legend.key.width = unit(1,"cm"), legend.margin = margin(rep(-3, 4))) +
labs(color = "Samples points", linetype = "Regression lines", shape = "", x = "Age (Years)") +
guides(color = guide_legend(override.aes = list(alpha = 1, size = 1.5, linetype = c(NA, NA)), order = 1), shape = guide_legend(override.aes = list(size = 1.5, alpha = 1), order = 2), linetype = guide_legend(order = 3)) +
ggtitle(bb)
p_list[[bb]] = p
}
p_f = p_list[["Feces"]] + p_list[["Mouth"]] + guide_area() + p_list[["Nose"]] + p_list[["Right_Arm"]] + p_list[["Right_Hand"]] + plot_layout(guides = "collect", nrow = 2) + plot_annotation(title = yy, tag_levels = 'A') & theme(plot.tag = element_text(size = 10, face = "bold"))
print(p_f)
#set_panel_size(p = p_f, g = patchworkGrob(p_f), file = paste0( "../output/figures/alpha_regression_", aa, "_", yy, ".pdf"), width = unit(3, "in"), height = unit(2.1, "in"))
}
}
